/*
 * EBFV1_assemblyMatrices.cpp
 *
 *  Created on: 13/09/2012
 *      Author: rogsoares
 */

#include "EBFV1_elliptic.h"

namespace PRS{
	double EBFV1_elliptic::assemblyMatrices(TMesh *pTMesh_data){
		double startt = MPI_Wtime();														// #START COUNTING TIME
		int nedges, n;
		int dim = pTMesh_data->getMeshDim();
		int ndom = pTMesh_data->getNumDomains();

		static bool key=true;
		if (key){
			allocateMatricesAndVectors(pTMesh_data);
			calculateAuxiliaryDataForMatrices_EFG(pTMesh_data);
			createAuxiliaryMappingForSolution(pTMesh_data);// STOP();
		}
		pPPData->calculateTotalMobility(pTMesh_data);
		for (int dom=0; dom<ndom; dom++){
			//  Matrices E,F,G are assembled edge-by-edge.
			pTMesh_data->edge_list[dom].getSize(nedges,n);
			for (int i=0; i<nedges; i++){
				EdgeData* edata = pTMesh_data->edge_list[dom].getValue(i,0);				// get next domain's edge
				VertexData* vdata_0 = edata->vertex_0;										// edge's vertex 0
				VertexData* vdata_1 = edata->vertex_1;										// edge's vertex 1

				vdata_0->countedges++;
				vdata_1->countedges++;
				for(int j=0;j<dim;j++){
					vdata_0->vel[j] = .0;
					vdata_1->vel[j] = .0;
				}

				int ID_0 = vdata_0->ID_PETScOrdering;										// vertex 0 ID for matrix assembly
				int ID_1 = vdata_1->ID_PETScOrdering;										// vertex 1 ID for matrix assembly
				double MobIJ = 0.5*(vdata_0->lambda_total + vdata_1->lambda_total);			// average total mobility between vertex 0 and vertex 1
				divergence_G(pMatVec->G,ID_0,ID_1,MobIJ,edata->aux_Matrix_G);
				divergence_E(pMatVec->E[dom],dim,ID_0,ID_1,MobIJ,edata);
			//	gradient_F_omega(pMatVec->F[dom],dim,ID_0,ID_1,edata);
			}
			// Matrix F needs boundary (external and inter-domains) faces contributions.
			gradient_F_gamma(pMatVec->F[dom],pTMesh_data);
//			PETSc_assemblyMatrix(pMatVec->F[dom]);
			PETSc_assemblyMatrix(pMatVec->E[dom]);
		}
		PETSc_assemblyMatrix(pMatVec->G);
		assembly_EFG_freeNodes(pTMesh_data);
		key=false;

//		printMatrixToFile(pMatVec->G,"Matrix_G.txt");
//		printMatrixToFile(pMatVec->E[0],"Matrix_E.txt");
//		printMatrixToFile(pMatVec->F[0],"Matrix_F.txt");
//		STOP();
		double endt = MPI_Wtime();															// #END COUNTING TIME
		return endt-startt;
	}

	void EBFV1_elliptic::calculateAuxiliaryDataForMatrices_EFG(TMesh *pTMesh_data){
		int nedges, n, i, j, pos;
		int dim = pTMesh_data->getMeshDim();
		int ndom = pTMesh_data->getNumDomains();
		const double I2D[4] = {1.0,.0,.0,1.0};
		const double I3D[9] = {1.0,.0,.0,.0,1.0,.0,.0,.0,1.0};

		for (int dom=0; dom<ndom; dom++){
			//  Matrices E,F,G are assembled edge-by-edge.
			int flag_dom = pTMesh_data->getdomFlags(dom);
			const double *K = pSimPar->getPermeability(flag_dom);
			pTMesh_data->edge_list[dom].getSize(nedges,n);
			for (int edge_th=0; edge_th<nedges; edge_th++){
				EdgeData* edata = pTMesh_data->edge_list[dom].getValue(edge_th,0);				// get next domain's edge
//				double sign = (edata->vertex_0->ID > edata->vertex_1->ID)?-1.:1.;
//				for (i=0; i<dim; i++){
//					edIJ[i] = sign*edata->versor[i];
//				}
//				for (i=0; i<dim; i++){
//					edata->versor[i] = sign*edata->versor[i];
//				}

				// auxiliary data for matrix G
				double KL[3] = {.0, .0, .0};
				pos = 0;
				if ( pSimPar->is_K_Isotropic() ){
					for (i=0; i<dim; i++){
						KL[i] = K[pos]*edata->versor[i];
						pos += dim+1;
					}
				}
				else{
					for (i=0; i<dim; i++){
						for (j=0; j<dim; j++){
							KL[i] += K[dim*i+j]*edata->versor[j];
						}
					}
				}
				edata->aux_Matrix_G = dot(edata->Cij,KL,dim)/(edata->length);//*edata->numRC);
				//printf("versor: %f %f  KL: %f %f lentgh: %f\taux: %f\n",edata->versor[0],edata->versor[1],KL[0],KL[1],edata->length,edata->aux_Matrix_G);

				// auxiliary data for matrix E
				// product between two vectors: Lij*Lij => matrix 3x3
				int k = 0;
				//const double *Lij = edata->versor;
				double matLij[dim*dim];
				for (i=0; i<dim; i++){
					for (j=0; j<dim; j++){
						matLij[k++] = edata->versor[i]*edata->versor[j];
					}
				}

				const double *Identity = (dim==2)?I2D:I3D;
				double matSubtrac_ILij[dim*dim];
				for (k=0; k<dim*dim; k++){
					matSubtrac_ILij[k] = -0.5*(Identity[k] - matLij[k]);
				}

				double EA[dim*dim];
				for (i=0; i<dim*dim; i++){
					EA[i] = .0;
				}
				k = 0;
				pos = 0;
				if ( pSimPar->is_K_Isotropic() ){
					int pos1 = 0;
					int pos2 = 0;
					for (i=0; i<dim; i++){
						for (j=0; j<dim; j++){
							EA[k++] = K[pos1]*matSubtrac_ILij[pos2+j];
						}
						pos1 += dim+1;
						pos2 += dim;
					}
				}else{
					for (i=0; i<dim; i++){
						for (j=0; j<dim; j++){
							for (k=0; k<dim; k++){
								EA[pos] += K[dim*i+k]*matSubtrac_ILij[dim*k+j];
							}
							pos++;
						}
					}
				}
				edata->aux_Matrix_E = new double[dim*dim];
				for (i=0; i<dim*dim; i++){
					edata->aux_Matrix_E[i] = EA[i];
				}
			}
		}
	}

	int EBFV1_elliptic::allocateMatricesAndVectors(TMesh* pTMesh_data){
		if (!allocateMatVec){
			pMatVec = new MatVec_struct;
			int np = pMData->getNum_GNodes();
			int numGF = pMData->getNum_GF_Nodes();
			int dim = pTMesh_data->getMeshDim();
			int ndom = pTMesh_data->getNumDomains();
			pMatVec->ndom = ndom;

			// RHS vector
			ierr = VecCreate(PETSC_COMM_WORLD,&pMatVec->RHS); CHKERRQ(ierr);
			ierr = VecSetSizes(pMatVec->RHS,PETSC_DECIDE,numGF); CHKERRQ(ierr);
			ierr = VecSetFromOptions(pMatVec->RHS); CHKERRQ(ierr);
			// Solution vector
			ierr = VecDuplicate(pMatVec->RHS,&pMatVec->solution);

			// Auxiliary vector z for "Matrix-free opertion"
			ierr = VecCreate(PETSC_COMM_WORLD,&pMatVec->z); CHKERRQ(ierr);
			ierr = VecSetSizes(pMatVec->z,PETSC_DECIDE,pMatVec->F_nrows); CHKERRQ(ierr);
			ierr = VecSetFromOptions(pMatVec->z); CHKERRQ(ierr);

			// Matrices G, E, F and Product EF for all mesh nodes
			ierr = MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,np,np,80,PETSC_NULL,80,PETSC_NULL,&pMatVec->G);CHKERRQ(ierr);
			ierr = MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,numGF,numGF,80,PETSC_NULL,80,PETSC_NULL,&pMatVec->G_freeNodes);CHKERRQ(ierr);
			ierr = MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,numGF,numGF,80,PETSC_NULL,80,PETSC_NULL,&pMatVec->EF_freeNodes);CHKERRQ(ierr);

			pMatVec->E = new Mat[ndom];
			pMatVec->F = new Mat[ndom];
			for (int dom=0; dom<ndom; dom++){
				ierr = MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,np,np*dim,100,PETSC_NULL,100,PETSC_NULL,&pMatVec->E[dom]);CHKERRQ(ierr);
				ierr = MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,np*dim,np,100,PETSC_NULL,100,PETSC_NULL,&pMatVec->F[dom]);CHKERRQ(ierr);
			}
			pMatVec->E_freeNodes = new Mat[ndom];
			pMatVec->F_freeNodes = new Mat[ndom];
			pMData->rowsToImport(pTMesh_data,pMatVec->nrows,pMatVec->rows);
			pMData->createVectorsForMatrixF(pMatVec->F[0]);
			for (int dom=0; dom<ndom; dom++){
				PetscInt F_nrows = pMData->get_F_nrows();
				const PetscInt* F_rows = pMData->get_F_rows_ptr();
				const PetscInt* F_cols = pMData->get_F_cols_ptr();
				PetscInt E_nrows = pMatVec->nrows;
				const PetscInt* E_rows = pMatVec->rows;
				const PetscInt* E_cols = pMData->get_pos_ptr();
				if (!F_rows || !F_cols || !E_rows || !E_cols){
					throw Exception(__LINE__,__FILE__,"Null pointer(s)!");
				}
				ierr = MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,E_nrows,np*dim,100,PETSC_NULL,100,PETSC_NULL,&pMatVec->E_freeNodes[dom]);CHKERRQ(ierr);
				ierr = MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,F_nrows,numGF,100,PETSC_NULL,100,PETSC_NULL,&pMatVec->F_freeNodes[dom]);CHKERRQ(ierr);
			}
			ierr = MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,np,np,80,PETSC_NULL,80,PETSC_NULL,&pMatVec->EF);CHKERRQ(ierr);
			allocateMatVec = true;
		}
		return 0;
	}

	int EBFV1_elliptic::assembly_EFG_freeNodes(TMesh *pTMesh_data){
		int ndom = pTMesh_data->getNumDomains();
		int nrows_F = pMData->get_F_nrows();
		const int* rows_F = pMData->get_F_rows_ptr();
		int ncols_F = pMData->getNum_GF_Nodes();
		const int* cols_F = pMData->get_F_cols_ptr();
		int np = pMData->getNum_GNodes();
		int dim = pTMesh_data->getMeshDim();
		int nrows_E = pMatVec->nrows;
		const int* rows_E = pMatVec->rows;
		int ncols_E = np*dim;
		const int* cols_E = pMData->get_pos_ptr();
		int nrows_G = pMatVec->nrows;
		const int* rows_G = pMatVec->rows;
		int ncols_G = pMData->getNum_GF_Nodes();
		const int* cols_G = pMData->get_idxFreecols_ptr();
		static MatReuse cll = MAT_INITIAL_MATRIX;
		static bool key = true;
		for (int dom=0; dom<ndom; dom++){
			if (key){
				ierr = MatGetSubMatrixRaw(pMatVec->F[dom],nrows_F,rows_F,ncols_F,cols_F,PETSC_DECIDE,cll,&pMatVec->F_freeNodes[dom]);CHKERRQ(ierr);
			}
			ierr = MatGetSubMatrixRaw(pMatVec->E[dom],nrows_E,rows_E,ncols_E,cols_E,PETSC_DECIDE,cll,&pMatVec->E_freeNodes[dom]);CHKERRQ(ierr);
//			printMatrixToFile(pMatVec->E_freeNodes[dom],"Matrix_E_free.txt");
//			printMatrixToFile(pMatVec->F_freeNodes[dom],"Matrix_F_free.txt");STOP();
		}
		ierr = MatGetSubMatrixRaw(pMatVec->G,nrows_G,rows_G,ncols_G,cols_G,PETSC_DECIDE,cll,&pMatVec->G_freeNodes);CHKERRQ(ierr);


		cll = MAT_REUSE_MATRIX;
		if (key){
			ierr = VecCreate(PETSC_COMM_WORLD,&pMatVec->z); CHKERRQ(ierr);
			ierr = VecSetSizes(pMatVec->z,PETSC_DECIDE,np*dim); CHKERRQ(ierr);
			ierr = VecSetFromOptions(pMatVec->z); CHKERRQ(ierr);
		}
		key = false;
		return 0;
	}

	int EBFV1_elliptic::PETSc_assemblyMatrix(Mat A){
		ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
		ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
		return 0;
	}

	int EBFV1_elliptic::createAuxiliaryMappingForSolution(TMesh *pTMesh_data){
		std::multimap<int,int> mapID_idxs;
		std::multimap<int,int>::iterator mapID_idxs__Iter;
		int nrows, ncols;
		int ndom = pTMesh_data->getNumDomains();

		// map node ID to its index (For multidomain meshes, ID over sub-domains will have two or more indices assigned. THat's why we use std::multimap)
		for (int dom=0; dom<ndom; dom++){
			pTMesh_data->vertices_list[dom].getSize(nrows,ncols);
			for (int j=0; j<nrows; j++){
				VertexData* vdata = pTMesh_data->vertices_list[dom].getValue(j,0);
				mapID_idxs.insert( pair<int,int>(vdata->ID,j) );
			}
		}
		// Let's transfer solution from pMatVec->solution to vertices_list.
		int numLocalIDs;
		int* pFreeNodesIDs;
		pMData->getRemoteIDs(numLocalIDs,&pFreeNodesIDs);
		freeNodesIndices.allocateMemory(numLocalIDs,ndom);
		for (int dom=0; dom<ndom; dom++){
			int k = 0;
			for(int j=0; j<numLocalIDs; j++){
				int ID = pMData->get_PETScToApp_Ordering(pFreeNodesIDs[j]+1);
				mapID_idxs__Iter = mapID_idxs.find(ID);
				if (mapID_idxs__Iter != mapID_idxs.end()){
					freeNodesIndices.setValue(k,dom,mapID_idxs__Iter->second);
					k++;
				}
			}
		}
		mapID_idxs.clear();
		return 0;
	}
}
